#include "AnaDef.h"
#include "XGLUtils.h"
#include "AliTRDdEdxUtils.h"
#include "TTreeStream.h"

TRandom3 gran(1);

Double_t means[]= {5,  6};
Double_t sigmas[]={1, 1.2};

const Double_t EPSILON = 1e-12;

void demo(Double_t &trq, Double_t &invl, Double_t &lgq, Double_t &mf, Double_t &inm, Double_t &jp, Double_t &lr, const Int_t kk)
{
  const Int_t nn = 100;//6;

  Double_t a0[nn];
  Double_t inva[nn];
  Double_t lga[nn];
  for(Int_t ii=0; ii<nn; ii++){
    a0[ii] = gran.Landau(means[kk],sigmas[kk]);

    if(a0[ii]<EPSILON){
      printf("a0[%d] = %15f; means[%d] = %15f; sigmas[%d] = %15f;\n", ii,  a0[ii], kk, means[kk], kk, sigmas[kk]); exit(1);
    }
    inva[ii] = 1/a0[ii];
    lga[ii] = TMath::Log(a0[ii]);
  }

  trq  =            XGLUtils::TruncatedMean(nn, a0,   0,   0.6);
  invl =         1./XGLUtils::TruncatedMean(nn, inva, 0.4, 1);
  lgq = TMath::Exp( XGLUtils::TruncatedMean(nn, lga,  0,   0.6) );
  //-----

  /*
  Double_t par[]={TMath::Mean(nn, a0), TMath::RMS(nn,a0)};
  XGLUtils::MaximumLikelihoodFit(nn, a0, XGLUtils::Landau, 2, par);
  //printf("MLFit: %15f\n", par[0]);
  mf=par[0];
  */

  /*
  //-------
  XGLUtils::MaximumLikelihoodFit(nn, inva, XGLUtils::InvLandau, 2, par);
  inm=par[0];
  */

  //-------
  Double_t tp0=1;
  Double_t tp1=1;
  for(Int_t ii=0; ii<nn; ii++){
    tp0*=TMath::Landau(a0[ii], means[0], sigmas[0], kTRUE);
    tp1*=TMath::Landau(a0[ii], means[1], sigmas[1], kTRUE);
  }

  jp = tp1/(tp0+tp1);

  lr = -999; //TMath::Log10(tp1)-TMath::Log10(tp0);
}

void jdet(Double_t &tpc, Double_t &trd, Double_t &co, Double_t &ws)
{
  const Double_t mean = 5;
  const Double_t lpc = 5*0.06;
  const Double_t lrd = 5*0.1;

  const Int_t nn = 100;

  TVectorD spc(nn), srd(nn);
  for(Int_t ii=0; ii<nn; ii++){
    const Double_t ipc = gran.Landau(mean, lpc);
    const Double_t ird = gran.Landau(mean, lrd);

    if(ipc<EPSILON || ird<EPSILON){
      printf("%d %15f %15f null;\n", ii,  ipc, ird); exit(1);
    }

    spc[ii] = 1/ipc;
    srd[ii] = 1/ird;
  }

  tpc = 1./XGLUtils::TruncatedMean(nn, spc.GetMatrixArray(), 0.4, 1);
  trd = 1./XGLUtils::TruncatedMean(nn, srd.GetMatrixArray(), 0.4, 1);

  TVectorD sc(2*nn);
  Int_t nc = -999; 
  co = AliTRDdEdxUtils::CombineddEdx(1, nc, &sc, 0x0, nn, &spc, 0x0, nn, &srd, 0x0)/4.99121505966755041e+00;

  const Double_t wpc = 6.60e-2;
  const Double_t wrd = 1.07e-1;
  ws = (tpc/wpc/wpc+trd/wrd/wrd)/(1/wpc/wpc+1/wrd/wrd)/5.00657367938121123e+00;

  tpc /= 5.01175380879650056e+00;
  trd /= 4.99320276031000443e+00;
}

void Loop(TTreeSRedirector *outStream, const Int_t ntrk, const Int_t kk)
{
  //mf: mlfit
  //jp: joint probability
  //lr: likelihood ratio
  for(Int_t itrk=0; itrk<ntrk; itrk++){
    if(itrk%(ntrk/100)==0)
      printf("Entry: %d / %d\n", itrk, ntrk);

    Double_t trq=-999, invl=-999, lgq=-999, mf=-999, inm=-999, jp=-999, lr=-999;
    demo(trq, invl, lgq, mf, inm, jp, lr, kk);

    Double_t tpc = -999, trd = -999, co = -999, ws=-999;
    jdet(tpc, trd, co, ws);

    (*outStream)<<"tree"<<
      "trq="<<trq<<
      "invl="<<invl<<
      "lgq="<<lgq<<
      //      "mf="<<mf<<
      //      "inm="<<inm<<
      "jp="<<jp<<
      //      "lr="<<lr<<
      "mean="<<means[kk]<<
      "sigma="<<sigmas[kk]<<

      "tpc="<<tpc<<
      "trd="<<trd<<
      "co="<<co<<
      "ws="<<ws<<
      "\n";
  }
}

int main(int argc, char *argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=2){
    printf("argc!=2\n");exit(1);
  }

  XGLUtils::SetFitPrintLevel(-1);

  TTreeSRedirector *outStream = new TTreeSRedirector("resdemo.root");

  const Int_t ntrk=atoi(argv[1]);//500000;
  Loop(outStream, ntrk, 0);
  Loop(outStream, ntrk, 1);

  delete outStream;

  return 0;
}

